Numerical computation of resonance poles in scattering theory 
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Abstract 



We present a possible way of computing resonance poles and modes in scatter- 



ing theory. Numerical examples are given for the scattering of electromagnetic 



waves by finite-size photonic crystals. 



Resonance poles are the main quantity of interest in scattering theory |T],|2|]. They are 
poles of the scattering matrix, and may be considered as generalized eigenvalues to which 
generalized eigenmodes are associated. These poles are complex ones, i.e. they correspond 
to complex values of the energy in the scattering theory of Schrodinger equation and to 
complex frequencies for the scattering theory of Maxwell system. 

When dealing for instance with finite size photonic crystals, i.e. periodic arrangements 
of scatterers that are finite in at least one direction of space 0, one cannot use Bloch waves 
theory to compute the electromagnetic properties of the structure and one has to retreat to 
the computation of the scattering matrix, that is the operator S (k) such that U = S (k) U l 
where U l is the harmonic incident field, with wavevector k, and U is the total field. The 
scattering matrix writes S (k) = I^ + T (k) where T (k) is the so-called scattering amplitude. 
Let us assume that there exists some pole k p of T in some neighborhood V of the complex 
plane, then locally the scattering amplitude writes T (k) = ^- + T (k) where P p is a residu 
operator and T is holomorphic in V. Operator P p is a finite rank operator and its range is 
precisely the kernel of T _1 (k p ). It is define in an abstract way as the Cauchy integral: 



where integration takes place on a loop oriented in the direct sense enclosing the only pole 
k p . Another way of defining the projection operator P p is to define it as the following limit 
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The point of this note is to show that the first abstract definition (|XJ) can be turned 
into a useful numerical tool for both the computation of the value of the pole and of the 



residu operator, and hence the generalized eigenmode, whereas the second is useless. From a 
numerical point of view, we of course only deal with finite rank operators and the scattering 
amplitude admits a representation as an operator on £ 2 (Z) , that is as a matrix, in the usual 
meaning, acting on double complex sequences @,[J. Once this representation is given, the 
residu operator can be computed provided that a region of the complex plane containing 
only one pole can be precised. This means that it suffices to know the value of the pole 
with a very poor precision to be able to compute the residu operator, which is not the case 
when the second definition @ is used: in that last case numerical instabilities necessarily 
occur as it uses the product of a singular matrix by a term tending to zero, which is a 
very bad numerical situation. From a practical point of view, one has to define a path 
7 : t G [0, 1] — > 7 (t) G C whose graph is a loop enclosing k p and to compute numerically 
the integral T (7 (£)) 7' (t) dt for which any reasonable numerical method works. However 
a precise computation of the pole is useful when one wishes to compute a map of the 
electromagnetic field of the pole, for in that case a particular basis such as Hankel-Fourier 
series are used, i.e. the field is expanded on the basis ( Hn (k p r) exp (in$) ) @||- A 

possible way is to use a Miiller like algorithm || and to compute a zero of the determinant 
of T -1 (k). However this matrix is generally badly conditionned and a better idea is to 
compute the smallest eigenvalue of T _1 (k). This works well in case of a finite size crystal, 
but this is not always the case: for instance, when modelizing photonic crystals by stacks of 
gratings 0] and introducing periodic defects, convergence problems may occur when using 
Miiller algorithm ||. We suggest then to compute the following Cauchy integral: 
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zT (z) dz = k p P p . (3) 



Recalling that P p has finite rank and hence has only a finite number of eigenvalues, a simple 
comparison of this last integral with P p gives the value of k p with a very good accuracy. 
Of course formula (f|) only holds when k p is a pole with multiplicity 1: this case is a very 
common one. One should not mistaken the range of P p and the multiplicity of k p : it is 
possible that the multiplicity of k p is 1 while the rank of P p is greater than 1 ||. 

Let us now turn to some numerical applications. We deal with the structure depicted 
in figure 1. It is a collection of 7 x 7 homogeneous fibers with relative permittivity e = 9, 
the radius of the rods R = 1/2 and the spacing is d — 1 (these values are given in arbitary 
units). We use a rigorous modal theory of diffraction to compute the scattering matrix of this 
system |4]|J. All the numerical results have been obtained using a standard PC computer. 
Removing a rod at the center of the crystal, a defect mode appears in the gap (see fig. 2 
for the transmission spectrum). To this peak in the transmission spectrum corresponds a pole 
k p . The reference value that we use for convergence comparison is k p = 2.32919703586134 — 
0.003782679876141 This value has been computed using Miiller algorithm by minimising the 
smallest eigenvalue of T _1 . For this given value of k p the smallest eigenvalue has modulus 
inferior to 10~ 14 . We then compute both Cauchy integrals (jl]j3|). We use an integration path 
that is a triangle whose vertices have affixes: (2.3,2.4 — 0.1z,2.4). We use the integration 



algorithm described in [11] and we denote by kw, P the numerical value obtained by using N 
points of integration, which we compare with the above value k p which is the best numerical 
value that we can obtain . A very good precision is rapidly obtained (see fig. 3): for 
instance, using a discretisation of 15 points we obtain 6 exact figures though with such a 
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rough discretization we only get operator P p with a low precision. In fact, it seems that 
the proportional coefficient between both integrals is not much affected by the precision 
with which P p is computed. A finer computation of integral (|1|) gives the defect mode. The 
convergence can be checked by looking at the non-zero eigenvalue of P p (fig. 4), here the 
reference eigenvalue (i.e. the best numerical value for a precision error of 10~ 15 ) is obtained 
with iV = 150. A much finer discretization than in the case of the pole is required to get a 
good representation of the defect mode, though the computing time is perfectly accessible 
with a very basic PC. 

In conclusion, we have shown that it was possible to turn a rather abstract mathematical 
object into a useful numerical tool. This technique applies as well for any situation in which 
a meromorphic operator with non essential poles is involved, which is the usual case. 
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Figures captions 

Figure 1: Sketch of the 2D photonic crystal. The transmission ratio is computed as the 
flux of the Poynting vector through the segment indicated below the crystal. 

Figure 2: Transmission ratio versus the wavenumber for an incident plane wave. 

Figure 3: Convergence of the value of the pole versus the number of integration points. 

Figure 4: Convergence of the eigenvalue of the projection operator versus the number of 
integration points. 
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